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We present results for the 1 = 2 tttt scattering length using Nf = 2 + 1 + 1 
twisted mass lattice QCD for three values of the lattice spacing and a range 
of pion mass values. Due to the use of Laplacian Heaviside smearing our 
statistical errors are reduced compared to previous lattice studies. A de¬ 
tailed investigation of systematic effects such as discretisation effects, vol¬ 
ume effects, and pollution of excited and thermal states is performed. After 
extrapolation to the physical point using chiral perturbation theory at NLO 
we obtain = -0.0442(2)stat(lo)sys- 
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1. Introduction 


Quantum Chromodynamics (QCD) describes, beyond the mass spectrum of stable and 
unstable hadrons, also their interaction due to the strong force. While originally thought 
to be not accessible to lattice QCD [T], Martin Liischer was able to relate the energy 
spectrum in finite volume to the infinite volume interaction properties in a series of 
seminal papers [a El n [5]. During the last years, extensions of these finite volume 
techniques have been worked out and they have been applied to a number of interesting 
systems. 

Liischer’s original formalism dealt with two massive bosons in a finite cubic box below 
the inelastic scattering threshold. It was originally formulated for the center of mass 
frame. In order to enlarge the applicability of Liischer’s method, extensions have been 
developed over the years which include: moving frames with non-vanishing center of mass 
momentum 1313 Eli, asymmetric boxes US [HIDE], twisted boundary conditions m 
naiiaiie], and more. Most of these aim to circumvent the poor resolution in momentum 
space due to the quantisation of the three-momentum in finite volume. Generalisations 
beyond the inelastic threshold for the two particle cases are also discussed [nuBiiini 
Eg ED |22] and various groups are now working on the more difficult three-particle 
scenario [231E1E51E61E3EE1E9]. 

In this first paper of a planned series of papers we investigate vr-Tr scattering in the 
isospin-2 channel. Trvr scattering in the 1 = 2 channel is technically the easiest case 
to consider in lattice QCD, because so-called fermionic disconnected contributions are 
absent. Moreover, the predicted quark mass dependence in chiral perturbation theory 
(yPT) at next-to-leading order (NLO) is governed by only one low energy constant, 
and so far lattice data showed surprisingly little deviations even from the leading order 
(LO) xPT predictions, which is parameter free. Isospin-2 mr scattering is, therefore, an 
important benchmark system to compare with non-lattice methods and an important 
test case for future investigations of hadron-hadron interactions, while simultaneously of 
phenomenological interest. 

Consequently, it has been computed in lattice QCD previously [30l ED ED ESI El] 
using Nf = 2 or Nf = 2-\-l dynamical fermions, for a recent review see Ref. j35|. In 
this paper we extend this list by using Nf = 2 -|- 1 -|- 1 dynamical quark flavours for 
the first time. The analysis relies on Wilson twisted mass fermions [36] at maximal 
twist [37] based on gauge configurations provided by the European Twisted Mass Col¬ 
laboration (ETMC) |38l 139] . This allows us to investigate a range of pion masses from 
250 to 500 MeV and discretisation effects by using three values of the lattice spacing. 
Eor estimating finite volume effects we have several ensembles at our disposal with all 
parameters fixed but the volume. 

We apply different strategies to determine scattering parameters from the data using 
Liischer’s finite volume method. This gives us an estimate of residual finite volume effects 
on the scattering length for other ensembles where we do not have multiple volumes at 
hand. 

The main result is an extrapolation of Mj^ao to the physical point utilising chiral 
perturbation theory at next-to-leading order. We carry systematic effects from different 
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ensemble 

/3 

a^ii 

ana 

afis 

{L/af X T/a 

A^conf 


^30.32 

1.90 

0.0030 

0.150 

0.190 

32^ X 64 

280 

^40.20 

1.90 

0.0040 

0.150 

0.190 

20^ X 48 

553 

A40.24 

1.90 

0.0040 

0.150 

0.190 

243 X 48 

404 

A40.32 

1.90 

0.0040 

0.150 

0.190 

32^ X 64 

250 

AGO.24 

1.90 

O.OOGO 

0.150 

0.190 

243 X 48 

314 

A80.24 

1.90 

0.0080 

0.150 

0.190 

243 X 48 

306 

A100.24 

1.90 

0.0100 

0.150 

0.190 

243 X 48 

312 


B35.32 

1.95 

0.0035 

0.135 

0.170 

32^ X 64 

250 

B55.32 

1.95 

0.0055 

0.135 

0.170 

32^ X 64 

311 

R85.24 

1.95 

0.0085 

0.135 

0.170 

32^ X 64 

296 

D45.32SC 

2.10 

0.0045 

0.0937 

0.1077 

32^ X 64 

301 


Table 1: The gauge ensembles used in this study. For the labelling of the ensembles we 
adopted the notation in Ref. [38]. In addition to the relevant input parameters 
we give the lattice volume and the number of evaluated configurations, A^conf- 


fit ranges through the entire analysis chain. In addition we investigate the stability of 
the extrapolation by different cuts in the pion mass values. 

The paper is structured as follows: in section we review the fermion action, the 
Laplacian-Heaviside smearing technique and the interpolating operators we use to con¬ 
struct the correlation functions. Liischer’s hnite volume method is explained in section]^ 
Our analysis and the treatment of all sources of systematic effects are worked out in sec¬ 
tion [4] The hnal results are discussed in section [5] and summarised in section [G] 

2. Lattice action 

The sea quarks are described by the Wilson twisted mass action with Nf = 2-|-l-|-l 
dynamical quark flavours. The Dirac operator for the light quark doublet reads |36j 

Di = D^^ + mo + , ( 1 ) 

where Dw denotes the standard Wilson Dirac operator and the bare light twisted 
mass parameter, and in general r*,z = 1,2,3 represent the Pauli matrices acting in 
flavour space. Di acts on a spinor xe = {u, and, hence, the u (d) quark has twisted 
mass {—fii)- 

For the heavy unitary doublet of c and s quarks m the Dirac operator is given by 

Dh = Dw + mo + ■ (2) 

The bare Wilson quark mass mo has been tuned to its critical value im EH). This 
guarantees automatic order O (a) improvement m, which is one of the main advantages 
of the Wilson twisted mass formulation of lattice QCD. 
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The splitting term in the heavy doublet Eq. introduces flavour mixing between 
strange and charm quarks which needs to be accounted for in the analysis. However, 
this is only important for quantities involving valence strange and charm quarks. 

2.1. Stochastic Laplacian Heavyside Smearing 

On the long term we plan to address as many observables as possible. Hence we have 
decided to use the stochastic Laplacian Heavyside smearing (sLapH), which is quite 
generally applicable and which is described in detail in Refs. |42l I43j . It represents a 
smearing method based on the covariant 3-dimensional Laplace operator 
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A“'(x, y; U) = (x)<5(y, x + k) + x - k) - 26{x, , (3) 

k=l 


where a, b are colour indices, x, y space-time coordinates, U denote smeared gauge links, 
and the sum is taken over the three spatial directions. The gauge fields in the Laplace 
operator, Eq. are smeared with three steps of three dimensional HYP smearing [U] 
with smearing parameters ai = 0:2 = 0.62, independently of the volume and lattice 
spacing. The Laplace operator can be decomposed as follows 

A = HaAaHI , (4) 

where Va represents the matrix of all eigenvectors and Aa is a diagonal matrix containing 
the eigenvalues. Colour, Dirac and space-time indices are suppressed. The smearing 
matrix then reads 

5 = ^ (5) 

where I 4 only contains eigenvectors corresponding to eigenvalues smaller than a cut-off 
al- These eigenvectors span what is called the LapH sub-space. The value can be 
chosen in a way that excited state contaminations are maximally suppressed. In Ref. [l3] 
it has been found that = 0.33 is optimal for excited state suppression independent of 
the interpolating operator. On a 48^ x 96 lattice, for example, this amounts to about 
960 eigenvectors per time slice. However, our tests show that it is better to have less 
eigenvectors for these large lattices, namely 660 on 48^ x 96 and 220 on 32^ x 64, to 
suppress excited state contaminations at early Euclidean times. 

The building blocks of observables are quark lines, Q, which can be written in the 
LapH framework as 

Q = sn-^S = H, , (6) 

where the middle part, V = is called perambulator with D = 74M and M 

being the Dirac operator. To generate a perambulator it is necessary to invert the Dirac 
operator for every eigenvector in 14 - It follows that for an all-to-all perambulator, since 
the Laplace operator is diagonal in time and Dirac space, this procedure needs to be 
repeated for each time slice and Dirac component. Taking the 48^ x 96 lattice as an 
example this would result in a total number of 253440 inversions. This number can 
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be reduced significantly by combining the LapH method with a stochastic approach as 
described in Ref. [13] in detail. 

For constructing all-to-all propagators with a stochastic approach, random vectors p 
are introduced which carry Dirac, time, colour and spatial indices. A propagator can be 
computed by solving 

^Xrlb] = prlb] ( 7 ) 


for Here the random vectors are not chosen to be random valued in all elements, 

but diluted, which is used to reduce the stochastic noise, see Ref. |l3| for details. Corre¬ 
spondingly, the composed index r[b] counts the total number of random vectors, Ar, via 
r and the total number of dilution vectors. Ad, via b. When constructing the all-to-all 
propagator via 




^ r=l b 


( 8 ) 


the zeros in the diluted random vectors ensure exact zeros in the product 
which reduces noise significantly. However, diluting random vectors leads to higher 
computational costs due to solving Eq. for each of these vectors separately. It is 
expected that the noise in correlation functions built from diluted stochastic propagators 
reduces with V\/iVR and ^/No, which favours more dilution vectors over more random 
vectors. Aside from this, each quark line needs its own set of random vectors to avoid a 
bias in the correlation functions. 

In the stochastic version of LapH, denoted by sLapH, the random vectors are intro¬ 
duced in time, Dirac and Laph sub-space indices. A quark line can then be estimated 
stochastically via 

E• (9) 

^ r=l b 


In our implementation we choose Z 2 random numbers. As mentioned before each quark 
line as defined in Eq. needs its own random vector to avoid a bias. This means that 
at least four random vectors are needed to be able to compute the correlation functions 
relevant for Trvr scattering processes. However, we use five random vectors, since the fifth 
random vector will allow for additional permutations in the four point functions. This 
improves the signal-to-noise ratio by a factor of 2 by only increasing the computing costs 
for inversions by 25%. 

As dilution scheme we have chosen full dilution in Dirac space combined with a block 
dilution in time and an interlace dilution in the LapH sub-space. The number of dilution 
vectors are summarised in table 2.1 They were chosen in such a way that the number 


of inversions and the noise in our observables are minimised simultaneously. 

We remark here that with the sLapH smearing scheme as explained above, only 
smeared-smeared correlation functions can be computed. Hence, we cannot compute 
matrix elements of local operators needed for instance for without major additional 
effort. 
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(L/a)^ X T/a 

iVD(time) 

A^d (Dirac) 

A^D(LapH) 

total ^ inversions 

20^ X 48 

24 

4 

6 

5 • 576 = 2880 

243 X 48 

24 

4 

6 

5 • 576 = 2880 

32^ X 64 

32 

4 

4 

5 • 512 = 2560 

48^ X 96 

32 

4 

4 

5 • 512 = 2560 


Table 2: Summary of the number of dilution vectors, A^D; used in each index. We use a 
block scheme in time and an interlace scheme in eigenvector space. The total 
number of inversions is the number of random vectors, here 5, multiplied by 
the number of inversions for one quark line. 

2.2. Operators 

The phase shifts are extracted from the finite-volume energies as described in the next 
section. In order to estimate the finite-volume energy spectrum, we build a matrix of 
correlators with a set of operators that resemble the vrvr isospin-2 system, and then use 
the variational method [13 B] to extract the spectral information. 

The SO{3) symmetry in continuum space is reduced to the octahedral group O on 
the lattice. For non-zero momentum P the symmetry is further reduced to the little 
group LG{P) that leaves P invariant. We construct the operators with definite total 
momentum P in each irreducible representation (irrep) of LG{P) via: 

C(.P,A, A;pi,p2)7r+(pi,t)7r+(p2,t + 1) • (10) 

Pl+p2=P 

pie{pi}* 

p2&{p2}* 

Here, 7r+(p, t) = t) 75 u(T, t), is the single pion operator projected onto mo¬ 

mentum p. With periodic boundary conditions, p is quantised as p = where n is 
a vector of integers. A is an irrep of the group LG{P) and A is the irrep row. {^ 1 , 2 }* 
represents the set of vectors {Rpi^ 2 -,R £ O}. C{P, A, X;pi,p 2 ) are the Clebsch-Gordan 
coefficient for Ai (8) A 2 —)• A, where Ai(A 2 ) is the irrep that 'n'~^{pi){T^~^{p 2 )) resides in. 
The two pions are separated by one time slice in order to avoid the complications due 
to Fierz rearrangement [46] . 

In this work we focus on building the operators for total zero momentum P = [0,0,0] 
and Ai(£ = 0,4,...) irrep where we can safely ignore all partial waves higher than 
^ = 0 in our analysis. The Clebsch-Gordan coefficients C(P, A, X;pi,p 2 ) are taken from 
Ref. [33] . 

The correlation matrix is computed via: 

cJ^(t) = (0|Op(t)O°’^(0)V)- (11) 

The variational basis contains various combinations of pi and p 2 that are allowed by the 
decomposition Ai (8) A 2 —)• A. For the case of P = [0, 0, 0] we include for the Ai irrep 
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the operators with |pi, 2 | = (0,1, 2, 3,4), which gives a 5-dimensional correlation matrix. 
More details about the irreps and other P values can be found in Ref. |34j . 

The energies are obtained by solving the generalised eigenvalue problem: 

C{t)vi{t, to) = Xi{t, to)C{to)vi{t, to). (12) 

It can be shown that the eigenvalues Xi{t) behave like 

+ , (13) 

where Ei is the i-th eigenvalue of the Hamiltonian of the system. However, we focus on 
zero total momentum where apparently solving the GEVP does not give any advantage 
over using pi = P 2 = 0 only. This is due to a rather weak coupling of different momenta 
in the matrix of correlators. All results presented in the following are, therefore, obtained 
directly from the four-point function at pi = P 2 = 0- 


2.3. Removing Thermal States 

As discussed in Ref. m and references therein, the spectral analysis of the two pion 
correlation function in the case of total zero momentum deviates from the usual cosh 
like behaviour. It was shown in Ref. [31] that the diagonal elements of the correlation 
matrix, Eq. [m for the two particle system obey the following spectral decomposition in 
the limit of large Euclidean times 

C^^(t) = Cg^^it) oc Aocosh(E„^(t - T/2)) + cexp(-M^T), (14) 


with Ej^j^ the two pion energy and the single pion mass, respectively, and constants 
Ao and c. This spectral decomposition differs from the standard by the term constant 
in Euclidean time. In the thermodynamic limit T —)• oo this polluting term vanishes, 
but for finite T it will dominate the correlation function at large Euclidean times. To 
remove this pollution, it was proposed to take a hnite difference first m and then build 
the following ratio m 


R{t + 1/2) 


Cwnit) — CT^T^it -\- 1 ) 

Cl{t)-Cl{t + l) 


(15) 


with C 7 r(t) the single pion two-point correlation function. One can show that the ratio 
has the functional form ini 


R{t + 1 / 2 ) = A(cosh(5E t') -|- sinh((I£' t') coth(2£'^f')) (16) 


with t' = t -\- II 2 — Tl2 and 5E = Etttt — 2M,r the energy shift. 

The generalisation of this procedure to correlation matrices and the variational method 
persists in shifting the correlation matrix Eq. before using the standard GEVP pro¬ 
cedure, see Ref. |3l] for details. 
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3. Finite Volume Methodology 


We are interested in the limit of small scattering momenta for the TTTr system with 1 = 2 
below inelastic threshold. Using the finite range expansion, the scattering length oq and 
the effective range rg can be related to the energy shift 5E by an expansion in 1/L as 
follows [3] 


6E 


dvroo 

Wj? 


+ Cl 


T 


+ C2^2 


+ C3 



Stt^Oq 


ro + 0(L-^), 


(17) 


with SE = i^TrTr ~ 2M^ and coefficients [a EH] 


Cl =-2.837297, C 2 = 6.375183, C 3 =-8.311951. 


More generally, including also non-zero total momentum, Liischer’s method relates the 
phase shifts 6 to the finite volume energy shift via the relation 


det 


-{M + i) 


= 0 , 


(18) 


where the matrix elements of the matrix M are given as [ 6 ] 

- 1 —iV jj 

— 1 ^ 3/2 ~j+i -2'js(l) Q )Clni,js,l' 


j=\l-l'\ s=-j 


(19) 


Z is Liischer’s generalised 2^-function and 


Q = 


27r 


( 20 ) 


the lattice scattering momentum. The elements of the tensor Cim,js,i'm' can be found 
in Ref. [ 6 |- For the case i = 0 with total zero momentum and no mixing with higher 


partial waves, Eq. 18 reduces to 


q cot 5o = 


L^/n 


ReizUhf)}- 


( 21 ) 


The scattering momentum q is then given as 

- Mj , (22) 

from which its lattice version, Eq. can be computed. 

As proposed in Ref. [ 6 | the continuum dispersion relation can be replaced by a lattice 
modified one. However, we do not see any difference in using one or the other version 
of the dispersion relation for the case of zero total momentum studied here. Hence, we 
stick to the continuum dispersion relation for this paper. 
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Figure 1: Ratios R(t) as defined in Eq. 15 as a function of tja for different ensembles in 


the center of mass frame. We also show the best fit curves for some represen¬ 
tatively chosen fit ranges. 


4. Results 

Before coming to the actual results, let us hrst describe our analysis procedure. Statis¬ 
tical errors are always computed using a bootstrap procedure with R = 1500 bootstrap 
samples. The bootstrap analysis is chained such that also statistical errors for best fit 
parameters can be determined. The gauge configurations are sufficiently separated HMC 
trajectories that autocorrelation does not play a role here, as we explicitly checked using 
the blocked bootstrap procedure. 

Energy shifts are determined from fully correlated fits of the ratios, Eq. [^to the data. 
The single pion energy levels needed for the fit are determined from the corresponding 
two point functions using a one state exponential (cosh) ht to the data, again fully 
correlated. The correlation between the single pion energy level and the ratio data is 
taken into account by using the same bootstrap samples. 
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In figure we show example plots of the ratio for various ensembles. For a represen¬ 
tative choice of the fit range we also show the corresponding best fit curves obtained 
by fitting Eq. 16 to the data. For the fit ranges shown one observes visually a good 
agreement between fitted curve and data. 

The fits to the ratios and the two point functions are repeated for a large number of 
fit ranges. We then assign a weight 


wx = {(1 - 2\px - 0.5|) • min(Ax)/Ax}^ 


(23) 


to every of these fits and quantities X = where px is the p-value of the corre¬ 

sponding fit and Ax the statistical error of {X) determined from the bootstrap procedure. 

Using the weights wx, we compute the weighted median over all fit results on the orig¬ 
inal data to obtain our estimate for the expectation value {X). The 68.54% confidence 
interval of the weighted distribution provides an (not necessarily symmetric) estimate 
for the systematic uncertainty stemming from the different fit ranges. The statistical 
error on {X) is computed from the R bootstrap samples of the weighted median. 

When derived quantities like q cot 5 are being determined, we follow the same proce¬ 
dure, just that the weights are now given by the products of the weights of the different 
contributing energy levels. 

The fit ranges are chosen such that for both the two point function and the ratio 
at least five time slices are included in the fit. The two point function is always fitted 
in an interval [tf,T/2] with > 6. For the ratio we fitted in the interval 
with E {11.5,13.5,15.5} and E {22.5,21.5,20.5,19.5,18.5} for L = 24,20 and 
E {30.5,29.5,28.5,26.5,24.5} for L = 32. We note that due to the weighting proce¬ 
dure described above we could have also included smaller values for q and tq without 
affecting the final result. 


4.1. Systematic Effects 

One of the main issues in this investigation is the question of systematic uncertainties 
in our analysis procedure. In particular, we have to consider possible contributions by 
the neutral pion which is much lighter than the charged pion on our ensembles due to 
twisted mass isospin breaking effects. Possible effects on the 1 = 2 pion scattering length 
extraction have been discussed in detail in Ref. |31j . 

In Ref. [3T] the authors were not able to find evidence of neutral pion contributions 
within their errors. Even though we have significantly smaller statistical errors, we 
still do not have any evidence for these effects within the statistical uncertainties. In 
the pseudo-scalar two point function the possible effect is an excited state with mass 
Mt^ + M^o due to the neutral pion having the vacuum quantum numbers in twisted mass 
lattice QCD [39]. This excited state could, however, never be identified. In the four- 
point function there can be effects from either close-by excited states at small Euclidean 
times or thermal states at large Euclidean times. Again, we do not see any evidence of 
those effects in our data. 

However, the analysis procedure detailed above is designed such that such effects 
should be covered by the systematic error we determine from the weighted distribution. 
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q; 
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o 
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m 

CD 


Histogram D45.32 


o 
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o 

lO 
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o 
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o 

CM 
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-0.312 -0.287 -0.262 -0.237 
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Figure 2: We show plain, weighted and bootstrap histograms for M-^ao and the ensembles 
A40.24 and D45.32, as explained in the text. In the lower left panel we show 
a QQ-plot for the bootstrap samples quantiles of M^-ao for A40.24. 
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ens 

a6E 

aq cot (5o 

AIt^CLq 

A30.32 

0.0037(1)(;^) 

-0.92(3)^^) 

-0.133(4)(+^) 

A40.32 

0.0033(l)(tJ) 

-0.90(3)e3) 

-0.155(5)(1^) 

A40.24 

0.0082(3)C^) 

-0.87(3)(tf) 

-0.164(5)e2) 

A40.20 

0.0179(5)(;i) 

-0.71(2)ej) 

-0.202(4)(+^) 

A60.24 

0.0076(2)(;j) 

-0.79(l)e|) 

-0.217(4)C2) 

A80.24 

0.0071(1)C?) 

-0.75(1)C0) 

-0.262(3)(1}) 

A100.24 

0.0063(1)C|) 

-0.75(1)C|) 

-0.294(3)(1?) 

B55.32 

0.0039(l)(t|) 

-0.71(2)ej) 

-0.219(5)C3) 

D45.32 

0.0084(2)(;°) 

-0.45(l)e°) 

-0.262(6)(;J2) 

B35.32 

0.0041(2)(;j) 

-0.82(3)e2) 

-0.151(6)e3) 

B85.24 

0.0085(1)C?) 

-0.66(1)C?) 

-0.292(3)(1^) 


Table 3: 6E, qcot5o and M-^qq computed with total zero momentum. 


This systematic error mirrors possible deviations from the theoretically expected curve 
and contributions by excited states or pollutions at large Euclidean times. 

The neutral pion might also contribute to the exponential finite size corrections in our 
data. As discussed below, for M,r and we take these twisted mass specific effects into 
account as determined in Ref. [50] from the data. For gcot((5o) finite size corrections 
specific to twisted mass have not been computed in yPT and, hence, we can only include 
the corrections computed in continuum yPT in Ref. m- However, these finite size 
corrections computed in continuum yPT have little influence on our results. Hence, we 
do not expect large effects from twisted mass specific hnite size effects. 

This conclusion is further supported by the fact that we do not observe large discreti¬ 
sation errors in the results for M^^aQ. Any contribution from the neutral pion should 
show up as a 0{a?) lattice artefact. Of course, all these indications are still not enough 
to finally exclude such systematic effects in our results. But we conclude that within our 
uncertainties they are negligible. 


4.2. Luscher formula to 0{1/L^) 


We will discnss several procedures to determine scattering parameters from the data. 


The first of which is to consider Eq. 17 to the order 1/L^ and, hence, neglecting the 
contribution from the effective range, oq can be determined from 5E, and L by 
numerically solving Eq. 17 for oq. The corresponding results for oq in units of M^- can 
be found in table [Sl 

As an illustration of our analysis procedure we show example histograms in fignre 
With plain histograms we mean histograms of the unweighted results for different fit 
ranges. In the weighted histograms the weights have been applied according to Eq. 23 


And with bootstrap data we denote the median of the weighted distribution evaluated 
on the bootstrap samples. In the histograms we plot the densities of the distribution. 
In the upper left panel the plain and the weighted histogram of M-^qq determined from 
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Figure 3: We show finite size corrected data for Mt^uq as a function of M -,^/with filled 
symbols compared to uncorrected data with open symbols. The corrected data 
include also the systematic uncertainty in order to allow for a comparison. 


Eq. 17 on the A40.24 ensemble are compared. The weighting leads to a well defined peak 
in the histogram, which is representative for the findings on most ensembles. 

In the upper right panel a comparison of the weighted histogram and the histogram of 
the bootstrap samples of Mt^clo again for A40.24 is shown. That the distribution of the 
bootstrap samples is approximately normal can be inferred from the lower left panel, 
where we show the QQ-plot of the bootstrap sample quantiles versus the theoretical 
quantiles of a standard normal distribution Here ft is the estimate of M^^ao 

and A/i its statistical error determined from the standard deviation over the bootstrap 
samples. This comparison indicates that the systematic and statistical uncertainties are 
approximately of the same size. This finding is again representative for most of the 
ensembles. 

In the lower right panel of figure we show again the plain versus the weighted 
histogram of Mn-oo; but this time for ensemble D45. This ensemble shows the largest 
systematic uncertainty of all the ensembles investigated here, which is also asymmetric. 
The weighted median and the systematic uncertainty which are indicated by the circle 
and the horizontal error bar above the histogram, respectively, show this asymmetry 
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clearly, even if it is not easy to identify by eye. 

For the successive analysis we also need to consider (exponentially suppressed) finite 
size corrections to our data. For and we use the results of Ref. [50] to correct 

our data. The corresponding data for f.,^ and the finite size correction factors are 
summarised in table in the appendix. We remark that A40.20 and D45.32 have not 
been considered in Ref. |50|. For D45.32 we use the factors computed in Ref. |50| for 
ensemble D20.48 with almost identical M 7 rL-value. For A40.20 we do not need the finite 
size corrections in the subsequent analyses. 

For ao we apply the asymptotic finite size correction formula Eq. (31) from Ref. [13| 

A(gcot(5) = {qcoi5)L - {qcot5)L=oo 

_-M^ ^ e—227 , \ (24) 

^ ^ 24nM^L + '7’ 

n=|n|7^0 

which is valid close to threshold (small q^) only. The c(n) are multiplicities for the 
ft and can be found for instance in Ref. m Note that at zero scattering momen¬ 
tum limq_^o 9 cot (5 = i/ao holds which gives us directly the finite size correction for the 
scattering length. 

For comparison we show the bare data for Mj^ao as a function of Mt ^/with open 
symbols and the corresponding finite-size corrected data with filled symbols in figure 
From this plot it is visible that most of the finite-size effects in this analysis stem from 
the ratio This is because the finite size corrections for M.^ and are opposite 

in direction. One can also see the effect of including the systematic uncertainties in 
the errors, which are included for the finite size corrected data points, but not for the 
uncorrected ones. This allows one to get an impression of their size and for which 
ensembles they are actually relevant. 

4.3. Liischer Formula to 0{1/L^) 

We have three A40 ensembles available, which differ only in their spatial extends, L = 
20, 24 and L = 32, respectively. Here, we can apply Eq. in order to estimate the 
scattering length oq and the effective range ro from the L-dependence. It amounts to a 
fit to the data for 6E for the three volumes with two fit parameters. Eor aM-,^ we used 
the values from the largest volume ensemble A40.32. The result is summarised in the 
first column of table [H 

According to the value (dof = 1) this is not a good fit. The data for 5E is shown 
together with the best fit and an error band in the left panel of figure]^ The result for 
is lower than the one for the A40.32 ensemble alone (see table [^, but it deviates 
not more than two a. The effective range parameter rg is determined by the fit only 
with large statistical uncertainties. 

It is likely that in particular A40.20 suffers still from exponential finite size artefacts. 
Therefore, we repeat the analysis with only A40.32 and A40.24 included. Now the fit 
basically reproduces the result to 0{1/L^) of the A40.32 ensemble, and there is no 
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1/L-fit 

1/L-fit 

q cot do-fit 

q cot do-fit 

L-values 
ao/a 

Al^do 

ro/a 

M^ro 

M^aoro 

xVdof 

32,24,20 

-0.98(5)(12g) 

-0.138(6)(12^) 

628(201) 

98(28) 

-87(23) 

5.14 

32,24 

-1.09(7)(120) 

-0.154(10)(1?8) 

42(221) 

6(31) 

-7(32) 

32,24,20 
-1.05(3)(l}i) 
-0.149(4)(1}6) 
147(22) (+30) 

21(3)(^^) 

-3.1(4)(;3) 

0.79 

32,24 

-1.09(6)(10) 

-0.154(8)(;J2) 

53(107) 

8(15) 

-1(2) 


Table 4: A summary of our extractions of the scattering length oq, and effective range 
ro for Liischer’s 1 /L expansion of the effective range expansion for two volume 
combinations each. 


sensitivity to the effective range parameter. The fit result is compiled in table and 
shown in the right panel of figure 


4.4. Finite Range Expansion for gcot^o 


Instead of using the expansion in 1/L as defined in Eq. 17, one may also determine 
q cot 6q close to threshold directly and perform the effective range expansion as follows: 


g'cotdo = — + , (25) 

ao 2 

where So is the s-wave phase shift, oq the corresponding scattering length, tq the finite 
range parameter and q the scattering momentum. Using the three ensembles A40.32, 
A40.24 and A40.20 again, we are able to obtain q cot So for three values of the squared 
scattering momentum q^. We then apply finite size corrections using Eq. 

The corresponding data are shown in the left panel of figure together with the best 
fit to expression Eq. The statistical error of the fit is indicated by the grey band. In 
addition, the value extrapolated to = 0 is shown. 

The best fit result is summarised in the third column of table [H In contrast to the fit 
of the 1/L expansion, the finite range expansion provides a good description of g cot do- 
The x^-value indicates a good fit, bearing in mind that there is little freedom left. The 
statistical uncertainties are smaller than the ones obtained from the 1/L expansion. In 
particular, a statistically significant value for the effective range parameter vq is obtained. 

Like in the previous case, where we used Liischer’s direct method, we also performed 
an extraction of uq and tq with only the largest two volumes L = 24,32 with the effective 
range formula. The result is summarised in the last column of table and the data are 
shown in the right panel of figure While the scattering length and effective range are 
in good agreement with previous extractions, we lose statistical significance for ro again. 
However, the error on oq is only slightly bigger than found for the other three fits. Note 
that the finite size corrections to q cot do do not change the fit result significantly. 
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a/L 


a/L 


(a) 


(b) 


Figure 4: (a) a5E as a function of a/L for the three A40 ensembles and the best fit 

(b) the same like (a) but using A40.32 and A40.24 only. 


according to Eq. 
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4.5. Chiral Extrapolation 

Because we have determined the scattering parameters for bare quark masses corre¬ 
sponding to larger than physical quark masses, we need to extrapolate to the physical 
point. For the / = 2 tttt scattering length yPT is particularly well suited, because it 
depends on only one low energy constant at NLO. 

As suggested in Refs. [521E3], it is convenient to write the yPT expression for 
as a function of because all quantities are dimensionless and no scale input is 

needed. The NLO yPT expression for written as a function of (Mjr//,r) at a xPT 

renormalisation scale = / 7 r,phys reads [52l [53] 




Ml f 

STTfl \ 167r2/2 


m2 

31n^ - 1 -inTrifJ-R 
J TT 


/TTjphys) 


(26) 


with related to the Gasser-Leutwyler coefficients as follows [5Tj 

Q '\ Gi 

innM = + -irh - h - 44 + 3 In . 

3 3 

Moreover, one can show in Twisted mass yPT that the leading lattice artefacts to 
are of order 0{a?Ml) [51]. At NLO we, hence, consistently describe our data with the 
continuum yPT formula provided above. 
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(a) (b) 

Figure 5: aq cot Sq as a function of for the three ensembles A40.32, A40.24 and 
A40.20 in (a) and for A40.32 and A40.24 only in (b). In addition the best fit 
of the effective range expansion to the data is shown as the line and the green 
triangle indicates the extrapolated value at = 0. The error band indicates 
the statistical nncertainty of the fit obtained from bootstrapping. 


Since we are not (yet) able to determine the pion decay constant /jr with the sLapH 
approach, we use data presented in Ref. m for the ratio //,r • The values for M^- //jr 
for all ensembles are compiled in table in the appendix. In this table we also give the 
valnes for the finite size correction factors as determined in Ref. m- The finite size 
corrections to oq are analytically compnted using Eq. 24 by setting g cot do = V“o- For 
the finite size corrections of in we also use the factors compiled in table as 

discussed above. 

The ratio can be determined with significantly smaller uncertainty than M-^aQ. 

Therefore, we do not expect that the missing statistical correlation between M .,^/and 
plays any role in the following analysis. 

For the fit we propagate the errors on and the finite size corrections on M-,^ 

and /tt nsing resampling. For the error on we add the statistical and systematic 

uncertainty of M^^ao and the statistical uncertainty of the finite size correction factor 
Kmt, in quadrature. In the minimisation we take errors both on M.,^ao and 
into account. We only include the large volume A40.32 ensemble (and not A40.24 and 
A40.20) in the fit. 

In the left panel of figure [^we show the fit to all the data (i.e. with no cut in M^//jr). 
The solid line represents the best NLO yPT fit to our data, the dashed line the LO 
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Figure 6: Left: Mj^oq as a function of M .,^/determined from Eq. 17 to 0{L with 


FS corrections for M .^/and Mj^aQ. The A-ensembles do not include A40.24 
and A40.20. We add the LO xPT prediction as the dashed line. The best fit 
to the data by the NLO xPT expression is shown as the solid line with error 
band. M^raolphys is plotted as the triangle. Right: the same but with LO xPT 
(Mjrao)'^® subtracted. 


parameter-free xPT prediction for M^oq. One observes that the LO xPT prediction 
already describes the data surprisingly well, and the NLO fit makes only a small cor¬ 
rection to the LO curve. We also show the value of M,rflo extrapolated to the physical 
point. In the right panel of figure we show the same, but with the LO xPT prediction 
(AfTrOo)^*^ subtracted. 

The data for the A- and B-ensembles fall to a good approximation on a single curve, 
whereas the only D-ensemble deviates slightly. While not inconsistent with a statistical 
fluctuation, this might be due to the rather small physical volume of the D45.32 ensemble. 
Another possible reason might be a?M^ lattice artefacts together with higher order terms 
in continuum xPT. 

We can explore this by restricting the fit range by applying an upper cut in the M,r //tt 
values. For the case of the fit including only data points with < 2.4 this is shown 

in figure which is otherwise identical to figure The cut value is indicated with 
brackets in the plots. 

The results of the fits to our data with different fit ranges are summarised in table 
We observe that different fit ranges do have little effect on the extrapolated value of 
Mjrflolphys- However, the fit quality improves with decreasing upper fit range, with the 
best fit for the cuts M,r// 7 r < 2.4 and M^/fn < 2.5. The only fit parameter shows a 
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p 

^'K'K 

■^Tr^^-olphys 


dof 

p-value 

M^Z/tt-cuI 

5.13(80) 

-0.0437(3) 

0.16 

2 

0.93 

2.2 

3.94(45) 

-0.0441(2) 

1.97 

4 

0.74 

2.4 

3.68(41) 

-0.0442(1) 

4.64 

5 

0.46 

2.5 

4.73(19) 

-0.0438(1) 

17.2 

8 

0.028 

- 


Table 5: Results from an NLO xPT fit to our data with different cuts in the upper 
fit-range. 




Figure 7: Like figure]^ but for the fit including only data points with M-n-ffTr < 2.4 as 
indicated by the brackets. 


large variation, though the fitted values are mostly consistent within errors. The reason 
for this large variation is visualised in the right panels in figures and [7| The curvature 
is mainly driven by the points around = 2. It is also visible that the fitted curves 

with cut still describe the data reasonably well within errors also for //jr values larger 
than the applied cut. 

As a final result we quote the weighted average of the two fits with cuts M^rZ/jr < 2.4 
and MttZ/tt <2.5, respectively 

M^ao = -0.0442(2)stat(^o)«ys> = 3.79(0.61)stat(lo:?l)sys, 

with the systematic error estimated from the maximal deviation of the single results in 
tableto the finally quoted one. The uncertainty stemming from the different fit-ranges 
of the fit to the ratio data turns out to be sub-leading. 
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5. Discussion 


The results presented in the previous section require some discussion. The first important 
point is that contributions from thermal and excited states as well as the contributions 
from the neutral pion can be excluded within the statistical uncertainty. With the use 
of the ratio, Eq. 16 thermal states constant in time are minimised if not cancelled 


completely. In section 4.2 e.g. figure we show histograms created from fits to a large 
number of fit ranges. Even for D45 where we see the largest systematic uncertainties 
stemming from the different fit ranges, the systematic uncertainty is not much bigger 
than the statistical uncertainty. This lets us conclude that we cannot resolve any of the 
aforementioned systematic uncertainties within our current statistical precision. 

Thanks to the three ensembles A40.32, A40.24 and A40.20 we were able to extract 
scattering parameters using three different methods: i) 1/L expansion, Eq. 17, to the 


order 1/L® for every A40 ensemble separately, ii) the same expansion but to the order 
1/L® and iii) a direct fit of the effective range expansion to the values of qcot6o{q‘^). 


Comparing the results in the sections 4.2, 4.3 and |4.4| shows that all three methods 


give compatible results for M-^ao, apart from method i) for A40.20. The reason for the 
deviation of method i) for A40.20 can be twofold: on the one hand = 2.99 might be 
too small and exponentially suppressed finite volume effects cannot be ignored. On the 
other hand the physical volume might be too small so higher orders in 1/L in Liischer’s 
formula and higher orders in in the effective range expansion are needed. Excluding 
A40.20 in methods ii) and iii) leads to smaller M-j^ao values, though still compatible 
within errors to the case when including A40.20. These smaller values, however, are in 
better agreement to method i) for both A40.32 and A40.24. From this fact alone it is hard 
to deduce which kind of volume effect we are dealing with here. That the physical volume 
of A40.20 might be too small is supported by the fact that the D45 ensemble, which is the 
other ensemble with small physical volume, gives a lower M-j^ao than expected although 
Mt^L = 3.87 is in the range of the other ensembles. The important consequence of this 
finding is that very likely all other ensembles do not suffer from significant finite volume 
effects. On every other ensemble than A40.20 and D45 the M^rL-value and the physical 
volume is at least as big as on A40.24 which seems to be free from volume effects within 
our current precision. 

For methods ii) and iii) also the effective range rg can be extracted from fits to the 
data. The best fit values vary significantly and have mostly very large uncertainties. 
Only for one fit (see table the ro-value is significant. We conclude that at the current 
level of available data and statistical precision we are not able to reliably extract the 
effective range parameter. 

The chiral extrapolation we present in section 4.5 includes three values of the lattice 
spacing. The A- and B-ensembles agree quite well, only at larger pion mass values small 
deviations show up. Unfortunately, we have currently only one ensemble for the smallest 
lattice spacing value available, namely D45.32. D45.32 has a quite small physical volume, 
while MttL ~ 3.8 which is compatible to the other ensembles. Currently, we cannot 
finally investigate the reason for this discrepancy. It could be a statistical fluctuation, 
which is supported by the exceptionally large systematic uncertainty for this ensemble. 
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Nf 

A/tj-G-o 

•^TTTT 

LO yPT 


-0.4438 


CGL (2001) 


-0.0444(10) 


CP-PACS (2004) 

2 

-0.0431(29)(-) 

- 

NPLQCD (2006) 

2+1 

-0.0426(6)(3) 

3.3(6)(3) 

NPLQCD (2008) 

2+1 

—0.04330(42)comb 

6.3(1.2)comb 

ETM (2010) 

2 

-0.04385(28)(38) 

4.65(0.85)(1.07) 

ETM (2015) 

2+1+1 

-0.0442(2) (1^) 

3.79(0.61)(lJ:?t) 

Yagi (2011) 

2 

-0.04410(69)(18) 

5.8(1.2)(-) 

Eu (2013) 

2+1 

-0.04430(25) (40) 

3.27(77)(-) 

PACS-CS (2014) 

2+1 

-0.04263(22)(41) 

- 


Table 6: Compilation of results for M^oo and including LO xPT, and 

Roy equations [56] denoted as CGL, CP-PACS |57|, NPLQCD (2006) [52|, 
NPLQCD (2008) |30j . ETM (2013) |31j . this work denoted as ETM (2015), 
Yagi et al. |58|, Eu |59| and PACS-CS [60] , 


or a finite volume effect, as discussed above. However, we can also not exclude a lattice 
artifact. We are working on additional D-ensembles, which will allow us to investigate 
this further. 

In table we provide a compilation of theoretical determinations for M^rOo and 
from the literature: this includes the LO xPT prediction as well as the value determined 
using yPT and Roy equations from Ref. [56| denoted as CGL. Eor the lattice results 
we have decided to include only direct determinations for which the chiral extrapolation 
has been performed: CP-PACS |57|, NPLQCD (2006) [52], NPLQCD (2008) [3^, ETM 
(2013) |3T|, this work denoted as ETM (2015), Yagi et al. [58], Eu [59] and PACS-CS [60] . 
We quote statistical and - where available - systematic uncertainties separately. Eor 
NPLQCD (2008) there is only the combined statistical and systematic uncertainty. 

The CP-PACS and the here presented ETM (2015) results are based on three values 
of the lattice spacing, two have been use by the authors of ETM (2013) and by Eu. 
The others have used only one value of the lattice spacing, however, NPLQCD (2008) 
and PACS-CS have used xPT to estimate discretisation effects. Einite size effects are 
estimated in the works of NPLQCD (2008), ETM (2013), Yagi et al. and Eu using chiral 
perturbation theory. In addition to the estimate from chiral perturbation theory we 
have studied in this work several volumes to obtain an estimate of residual finite volume 
effects. 

As can be observed visually in figure the values for show an overall good 

agreement. In particular, there is no definite dependence on Nj. The exception are 
the high values of NPLQCD (2006) and PACS-CS. PACS-CS quotes in addition a small 
uncertainty, which makes the deviation to our result statistically significant. If the 
systematic uncertainty is added to the statistical error, the results are compatible again. 
PACS-CS has studied the smallest pion mass value around 170 MeV of all the works 
considered in this comparison. However, they claim they cannot fit this point with NLO 
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Figure 8: Comparison of the various predictions for Mj^aolphys- The statistical error is 
indicated by the error bars on the points, the systematic uncertainty, where 
available, as coloured bars. Note that for NPLQCD (2008) statistical and 
systematic errors are combined. 


yPT. Therefore, they had to include higher order terms. More results at close to physical 
pion mass values might, therefore, clarify whether this is a statistical or systematic effect 
in the data, or whether the chiral extrapolation with NLO yPT is misleading. 

The values for show a large variation, which is, however, covered by the large 
errors on this LEG. 


6. Summary 

In this paper, low-energy tttt scattering in the isospin 1 = 2 channel is studied within 
Liischer’s finite volume formalism in lattice QCD. We use for the first time Nf = 2-|-l-|-l 
dynamical quark flavours based on gauge field configurations provided by ETMC. The 
list of ensembles covers three values of the lattice spacing, several volumes and a large 
range of pion mass values, see table We determine energy shifts 5E = — 2Mt^ 

using the stochastic LapH method and convert them in scattering length values applying 
Liischer’s formalism. 

We apply different, though closely related, methods to determine scattering param¬ 
eters and find compatible results for M^^aQ. The finite range parameter rg cannot be 
determined with sufficient certainty. 

Due to the use of ensembles with a broad range of pion masses we are able to extrapo¬ 
late Mj^ao towards the physical pion mass point. After correcting for finite volume errors 
of the scattering length oq we are able to make a rather smooth extrapolation towards 
the physical pion mass value, the result of which is shown in figure]^ We do not observe 
significant discretisation effects between A and B-ensembles. The only D-ensemble D45 
shows a deviation which can be equally well explained by discretisation effects, finite 
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volume effects or by a statistical fluctuation. To clarify this point, we are generating 
data on more D ensembles with smaller pion masses. 

Our final result for is, = —0.0442(2)stat(lo)sys- We have compared our 

result to a list of other lattice determinations and combined with Roy equations 
and found mostly remarkably good agreement. 

The currently ongoing extension of the study presented in this paper is mr scattering 
with 1=1, where the p resonance is present. With the perambulators ready from 
the distillation process, this is straightforward to do and first results are available. For 
the pion-pion scattering, the channel / = 0 is more complicated, particularly for the 
twisted mass formulation due to isospin breaking at finite lattice spacing values. Another 
currently ongoing extension is to address low-energy scattering of other mesons: e.g. ttK 
and KK scattering and more than two mesons. Such techniques can also be extended 
to charmed meson scattering processes relevant for the recently discovered XYZ states. 
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A. Data Table 




M^/U 



A30.32 

0.1239(2)(Zi) 

1.915(10) 

1.0081(52) 

0.9757(61) 

A40.32 

0.1415(2)(Z|) 

2.068(08) 

1.0039(28) 

0.9874(24) 

A40.24 

0.1446(3)(Z}) 

2.202(13) 

1.0206(95) 

0.9406(84) 

A40.20 

0.1474(6)(Zi) 

NA 

NA 

NA 

A60.24 

0.1733(3)(Z^) 

2.396(11) 

1.0099(49) 

0.9716(37) 

A80.24 

0.1993(2)(;i) 

2.623(07) 

1.0057(29) 

0.9839(22) 

A100.24 

0.2224:{2){tl) 

2.788(07) 

1.0037(19) 

0.9900(15) 

B35.32 

0.1249(2)(Z|) 

2.047(11) 

1.0069(32) 

0.9794(27) 

B55.32 

0.1555(2)(Z|) 

2.352(07) 

1.0027(14) 

0.9920(10) 

B85.24 

0.1933(3)(;?) 

2.736(15) 

1.0083(28) 

0.9795(24) 

D45.32 

0.1207(3)(;|) 

2.485(12) 

1.0047(14) 

0.9860(13) 


Table 7: Single pion energy levels, and the finite size correction factors and 

Kf^ computed in Ref. [50] for M.^ and /^, respectively. 


27 








